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ABSTRACT 

We present the results of numerical simulations of shock wave-driven jets in the solar atmosphere. 
The dependence of observable quantities like maximum velocity and deceleration on parameters such as 
the period and amplitude of initial disturbances and the inclination of the magnetic field is investigated. 
Our simulations show excellent agreement with observations, and shed new light on the correlation 
between velocity and deceleration and on the regional differences found in observations. 
Subject headings: magnetic fields — MHD — Sun: chromosphere — shock waves 
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1. INTRODUCTION 

In the solar chromosphere, jet-like dynamic features 
are found in several regions. At the quiet sun limb, 
we find spicules, chromospheric protrusions that reach 
heights of 5-9 Mm and la st for 3-15 mi nutes, reaching 
velocities of 10-30 km s _1 (|Beckerslll968[ ). Man y models 
have been proposed to explain their formation ([Sterling! 
2000), but observational and interpretational difficulties 
have made the models hard to constrain. 

On the quiet sun disk, we observe dark mot- 
tles, which appear to have many similarities to 
spicules. There has been some controversy over 
the relationship between spicules a nd mottles, with 
iGrossmann-Doerth fc Schmidtl (|l992l) concluding that 
they are not counterparts, while other authors have ar- 
gued that the similarities a re striking dTsiropoula et al.l 
119941 : ISuematsu etallll995t IChristopoulou et alJl2QQlf l. 

A third group w hich has c ome under study rec ently 
(|De Pontieu et al.ll2004l . 120071 : lHansteen et~alfj2006[ ) are 
dynamic fibrils, a subset of jets frequently found in the 
vicinity of active region plage. These jets are shorter than 
spicules in both length and duration, reaching heights 
of 1-4 Mm and lasting 3-6 minutes, and frequently 
show both periodicity and internal structure varying on 
shorter timescales. 

Recent observationa l evidence (Ha nsteen et al.l [2006; 
iDe Pontieu et al1l2007l ) has suggested that dynamic fib- 
rils are driven by magnetoacoustic shocks, which orig- 
inate in the convection zone and photosphere and, al- 
though usually evanescent in the chromosphere, may be 
able to leak into the upper layers in inclined or heated 
flux t ubes. The data of iRouppe van der Voort et al.l 
(|2007f ) indicate that the same mechanism may also be 
the driving force behind at least a subset of quiet sun 
mottles. 



We perform numerical simulations of such 
shock wave-driven jets (henceforth called fib- 
rils), and ana ly se th e data in a simila r fashi on to 
Hansteen et al.l (120061). IDe Pontieu et all (|2QQ7f ) and 
Rouppe van der Voort et aTl ( 20071 ). With a more 



idealised and controlled environment, we can study the 
effect of parameters such as the period and amplitude 
of the piston driver and the inclination of the magnetic 
field. We examine the correlation between the maximum 
velocity of a fibril and its decele ration, and also see 
if we can confirm t he suggest i on (|M ichalitsanos 1973; 
Bel fc Lerovl 119771 : Jguematsul fl990MDe Pontieu et al.l 
2004 120051 . 120071 : lHansteen etafl l2QQ6h that inclined 
magnetic fields can allow normally evanescent long- 
period waves to propagate into the upper chromosphere 
and coron a. Signs of such waves have recently been 
observed (IDe Pontieu et all 120051 : iMcIntosh fc Jefferiesl 
120061 : Jefferies et al. 20Qffy. 
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2. SIMULATIONS 

For our simulations, we use a simple one-dimensional 
model of the upper solar atmosphere, with a monochro- 
matic piston driver at the lower chromospheric boundary 
for creating acoustic waves. As these waves travel up- 
wards, they gain in amplitude because of the decreasing 
density of the medium and steepen into shocks (assuming 
a large enough initial amplitude), which hit the transi- 
tion region and thereby push the corona upwards. 

Fig. [T] shows the initial density and temperature pro- 
files of our model atmosphere. The model extends 
around 8.5 Mm in height, starting 0.85 Mm above the 
photosphere, with the transition region about 1 Mm 
above the lower boundary, the chromosphere below and 
the corona above it. The temperature at the upper 
boundary is maintained at 1 MK while heat conduction 
and radiation set the temperature in the rest of the do- 
main. 

In computing the radiative losses we include a ra- 
diative loss term that accounts for the energy loss due 
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Fig. 1. — Initial density (dashed, left axis) and temperature 
(solid, right axis) structure of the model. 

to collisional excitation of the various ions comprising 
the plasma. We have included the elements hydro- 
gen, carbon, oxygen, neon, and iron, as well as ther- 
mal bremsstrahlung, using the ionization and reco mbi- 
nation rates given by lArnaud fc Rothenflu j (p. 985) and 
Shull fc van Steenber and using the collisional 

excitation rate s found through the HA O-DIAPER atom 
data package ([Judge fc Meisneri [l994). The metals are 
treated by assuming ionization equilibrium and then de- 
riving an a priori radiative loss curve as a function of 
electron temperature. While one should ideally solve the 
equation of radiative transport in order to calculate the 
radiative losses from hydrogen, comparison with models 
where this has been done indicate that the errors in- 
curred by assuming effectively thin losses in the Lya line 
are not significant to the processes studied in this pa- 
per. In order to avoid too large radiative losses in the 
lower atmosphere the radiative loss term is multiplied 
with exp(— r), where the "optical depth" r is set propor- 
tional to the gas pressure. 

Thermal conduction is given by K T 5 / 2 dT/ds, which 
carries energy from the corona to the upper chromo- 
sphere and determines the temperature structure of the 
transition region. The numerical model is set up so that 
conduction is carried along the magnetic field. The en- 
ergy equation is solved by operator splitting with thermal 
conduction solv ed implici t ly u sing a multi-grid method 
as described by iHansteenl (|2QQ5f ). 

Although it might have been ideal to include the so- 
lar atmosphere from the photosphere to the corona in 
our model, we have chosen to restrict ourselves to the 
magnetically dominated upper regions, from the chromo- 
sphere and upwards. This is because, with an inclined 
magnetic field, the (3 « 1 layer, where (5 is the ratio 
between the thermal and magnetic pressure, is an area 
of extensive conversion between the so called fa st and 
slow modes (jHegglandl l200l iBogdan et aTll2003[ ). We 
are here primarily interested in the field-aligned acoustic 
(upper-atmosphere slow mode) waves, and not the mag- 
netosonic fast mode waves that can propagate across field 
lines. The reason is that, in a real multidimensional at- 
mosphere, the fast mode wa ves will be refract ed into ar- 
eas of low Alfven speed (e.g. IOsterbrocklll96ll ) , and thus 
have difficulty propagating into the higher layers of the 
atmosphere. In a ID model, they are forced to propa- 



gate vertically, and will then reach the upper layers with 
far greater energy than is realistic. For this reason, we 
get cleaner and arguably more realistic results by plac- 
ing our driving piston in the low-/? region of the (upper) 
chromosphere, and restricting our domain to the areas 
above that. 

Of course, (3 could be reduced to any value we like 
simply by increasing the strength of the magnetic field. 
We have chosen a field strength of 6.0 x 10 -3 T (60 
gauss), representing a reasonable value for chromospheric 
field outside of, but not too far from, active regions. A 
stronger field would allow us to include deeper layers of 
the atmosphere, but at the cost of shorter timesteps be- 
cause of the increased Alfven speed. We believe that 
including deeper layers would be unlikely to change our 
results significantly, especially in light of the other sim- 
plifications inherent in a ID model. In ID, we are pri- 
marily interested in studying fundamental effects, and 
further refinements should rather be considered in future 
multi-dimensional investigations. 

The code w e use is a version of the code described 
in IHansteenl (|2005l ). (See also Galsgaard & Nord- 
lund's description of an earlier version of the code at 
http : //www. astro .ku.dk/~kg). This code allows us to 
use 3D variables in a ID domain, thus allowing for trans- 
verse components of the velocity and magnetic field. We 
exploit this by using the same background atmosphere 
(density /temperature) in all cases but tilting the mag- 
netic field away from the vertical. 

With plasma movement thereby mostly constrained to 
being along the inclined magnetic field, but having a ver- 
tical computational domain, the simulated data are ba- 
sically projections of the field-aligned motion onto the 
vertical. As a result, propagation speeds appear lower 
when the field is highly inclined — at 60°, the waves 
have to travel twice as far along the field to reach the 
same height. If desired, the vertical velocities can be re- 
calculated as field-aligned by multiplication with a simple 
1/cos factor, being the inclination. 

Using a ID model makes our analysis simpler, while 
keeping the important physics. The main phenomena 
we exclude are wave refraction, which (as mentioned) is 
only important for the fast mode, and the curvature and 
expansion of the magnetic field. 

In our simulations, the 60 gauss field is constant and 
homogeneous. We vary the inclination, from 0° (verti- 
cal) via 30° and 45° to 60°. We use piston periods of 180, 
240, 300 and 360 s, and initial amplitudes of 200, 500, 
800 and 1100 m s _1 at the lower boundary. The driv- 
ing is sinusoidal and the piston is active throughout the 
simulations. The piston movement is along the magnetic 
field rather than vertical, in order to further suppress 
unwanted fast modes. 

The piston generates a train of waves that quickly 
(especially at greater initial amplitudes) steepen into 
shocks. As these shocks hit the transition region, they 
give a large impulsive acceleration (a "kick") to the 
plasma there as the shock front passes, then a more 
gentle "push" for a while in the receding phase of the 
shock. Some of the wave energy passes through and en- 
ter s_Jh^_ccTOna 2 _but quite a lot is reflected down again 
— iHegglandl (|2QQ3h estimates a reflection coefficient of 
about 0.70 for low frequency linear waves, though shock 
waves may behave slightly differently. 
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In Fig. [21 we have plotted a representation of the ver- 
tical energy flux for an example case as a function of 
height and time. Reflection from the transition region 
and re-reflection from the closed lower boundary lead to 
the formation of standing waves in the chromosphere, 
which makes reliable estimation of the upward energy 
flux difficult; however, crude estimates comparing the 
input energy from the piston with the energy flux reach- 
ing the upper boundary give a transmission coefficient of 
10% or less, the rest being reflected or dissipated. 

A plot of the simulated vertical velocity as a function of 
height and time for an example case is found in Fig. 03 In 
this case we hardly see the standing waves because of the 
large amplitudes the shock waves reach when they enter 
the transition region where the density falls off rapidly 
with height. 

3. ANALYSIS 

I n their observational st udy, lHansteen et al.l (2006) 
and lDe Pontieu et al.l (|2007l ) find that fibrils, which they 
observe in Ha, follow parabolic paths in distance-time 
diagrams. These parabolic paths are found in our sim- 
ulations as well, cf. Fig. [3j They can be vis ualized in 
a way more comparable to iDe Pontieu et all if we in- 
stead show the temperature in a distance-time plot, as in 
Fig. [H We clearly see the parabolic paths traced by the 
transition region as it is periodically pushed up by the 
passing shock waves. Ha radiation is primarily generated 
in the hot upper layers of the chromosphere just below 
the transition region, and observations in that band will 
therefore show very similar movement. 

A look at Fig. [3] tells us that, although the parabolic 
shapes are clearly seen, there are also other features 
present. Some fast mode waves remain — one can be 
seen faintly at the beginning of the simulation, before 
the first slow shock front — and some additional slow 
mode waves are apparently generated in the shock fronts 
through nonlinear processes. Although these other waves 
can have significant amplitudes, their effect on the move- 
ment of the transition region is far less than that of the 
slow mode shock waves. In addition, fast modes would 
tend to refract downwards in a real atmosphere, as dis- 
cussed in the previous section. 

The picture is clearer when looking at the temperature 
plot. Therefore, it seems reasonable to use the movement 
of the transition region as a proxy for the movement of 
the layers where the main Ha radiation takes place, make 
parabolic fits of the position, and use derivatives of these 
to determine velocities an d decelerations. This is simi- 
lar to the method used bv lDe Pontieu et al.l (|2007l ). and 
makes comparisons easier. Furthermore, it means that 
the calculated decelerations will be constants. An exam- 
ple of the parabolic fits is shown in Fig. [5l 

Having determined the decelerations and velocities of 
each fibril from our 64 simulated cases (not all of which 
lead to noticeable movement of the transition region), 
a total of 19 fibrils, we can make scatterplots in the 
same way as IDe Pontieu et al.l The results are shown in 
Figs. [6] and [3 — the latter is corrected for projection ef- 
fects through multiplication by 1/cos when the field 
is inclined. We see obvious linear correlations between 
deceleration and maximum velocity, but the slopes are 
clearly different in the 180 s and 240 s cases. Although 
the number of points for the 300 s and 360 s cases is 



2000 



1500 



V) 

qj 1000 
E 



500 





1.0 1.5 2.0 2.5 5.0 5.5 
Height [Mm] 

Fig. 2. — The vertical energy flux density of the waves in an 
example case (45° inclination, 240 s driver period and 1100 m s _1 
initial amplitude). We have plotted the square root of the actual 
energy flux density in order to preserve the sign of the velocity — 
white represents upward velocity and black downward. Most of the 
energy ends up in standing waves in the chromosphere, while only 
a small fraction (almost invisible in this plot) reaches the corona. 
The dotted line is at the height where the temperature is 100 000 K, 
i.e. the lower transition region. 
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Fig. 3. — The vertical velocity taken directly from the simulation 
data in the same case as Fig. [2] In addition to the main fronts and 
the parabolic shapes, we see an initial (quite weak) fast mode front 
and some additional fronts that seem to be generated at the tops 
of the fibrils. Again, the dotted line marks where T=100 000 K. 
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Fig. 4. — The gas temperature as a function of time and vertical 
position in the same case as for Fig. [3] This way of plotting very 
clearly shows the movement of the transition region which is our 
primary interest. 
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Fig. 5. — An example of the parabolic fitting process, for the 
same case as Figs. [3] and 3] The transition region height is defined 
as the height where the temperature is 100 000 K. 

small, they too appear to show the same trend, giv- 
ing progressively lower decelerati ons for a gi ven maxi- 
mum velocity. The data of iDe Pontieu et al.1 are much 
more scattered, both because of the greater difficulty in 
analysing multi-dimensional observational data and be- 
cause of uncertainty in the relative orientations of the 
fibrils' movement and the line of sight, but they also find 
a more or less linear correlation and signs of a period- 
dependent slope — again giving smaller decelerations at 
greater periods. 

A question of interest is how big a role gravity plays in 
setting the deceleration. One natural instinct when see- 
ing the parabolic paths is to assume that we are observ- 
ing a ballistic flight under the influence of solar gravity, 
with differences being accounted for by inclination and 
line of sight effects. In observations, the angle between 
the fibril movement and the line of sight is often difficult 
to determine exactly, but the many low decelerations re- 
ported would require extremely steep inclinations for the 
fibrils (jSuematsu et al.lll995l : foe Pontieu et al.ll2QQ7f ). 

In looking at our simulation data, we have no line of 
sight uncertainty, and compensating for the inclination 
of the fibrils is trivial. We can therefore easily calcu- 
late the distribution of decelerations with inclination an- 
gle. This distribution is shown in Fig. El The black 
histogram shows the field-aligned component of the so- 
lar gravitational acceleration, while the data points are 
shown with the same shapes and colour scheme as used 
in Fig. [U We see that all our observed decelerations 
are smaller than the projected gravity, and often signif- 
icantly smaller. Furthermore, they are not particularly 
clustered, pretty much the whole range of decelerations 
from about 50 m s _1 to gravity being represented at all 
angles. With all projection effects already accounted for, 
this distribution is extremely hard to reconcile with a 
ballistic model, and indicates that we should look not to 
gravity but elsewhere in trying to explain it. 

One alternative is to look at shock wave physics and 
pressure forces. The period of a wave train should re- 
main constant during its propagation, and in that case, 
a fully developed shock wave (an N-wave) has a given 
time during which the amplitude must move from its 
highest to its lowest value. If the maximum amplitude 
— approximately the same magnitude as the maximum 
velocity of the fibril — increases, the deceleration must 
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Fig. 6. — Scatter plot of the maximum (upward or downward) 
velocities of fibrils vs. their decelerations, showing clear linear 
correlations. Red triangles correspond to a period of 180 s, blue 
squares to 240 s, green crosses to 300 s and black diamonds to 
360 s. The slopes are clearly different as the period increases. 
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Fig. 7. — As Fig. [6] but adjusted for projection effects. 
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Fig. 8. — Distribution of decelerations with magnetic field incli- 
nation. The data symbols are the same as in Fig. [6] All deceler- 
ations are lower than the projected gravity (solid columns), often 
very much lower. Since we have already compensated for projec- 
tion effects, this distribution is extremely hard to reconcile with a 
ballistic theory. 
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Fig. 9. — An N-wave with lower amplitude (dotted line) will 
have a less steep descending slope (lower deceleration) than one 
with higher amplitude (solid line) at a given period. Similarly, one 
with a longer period (dashed line) will have a lower deceleration 
than one with a shorter period at a given amplitude. 
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Fig. 10. — As Fig. [3 but including the theoretical values from the 
shock deceleration hypothesis as solid lines. There is a near-perfect 
fit with the data, giving strong support to the hypothesis. 

then be greater per unit time. Similarly, if the period is 
longer, the wave has more time between maximum and 
minimum amplitude, and the slope becomes less steep, 
meaning the deceleration is lower for a given amplitude. 
This is illustrated in Fig. [9l 

The shock deceleration hypothesis can thus explain 
both the linear correlation between deceleration and 
maximum velocity and the variation in its slope with 
wave period. Indeed, we can calculate the decelerations 
expected from the theory as a function of period and 
maximum velocity /amplitude and compare them with 
the values from our simulations. The formula is simply 



d = 



P/2 ' 



(1) 



where d is the deceleration, v max is the maximum veloc- 
ity and P the period of the wave. As we see in Fig. [10l 
there is near-perfect correspondence between the theoret- 
ical and simulated values, giving strong support to this 
explanation. 

Although all the fibrils we study in this paper have de- 
celerations below projected gravity, the shock wave de- 
celeration model can also explain decelerations greater 
than gravity. Indeed, in preliminary experiments with 
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Fig. 11. — Temperature plots showing the movement of the tran- 
sition region with time. The period is 300 s; the piston amplitude 
is kept constant at 500 m s — 1 , while the inclination of the mag- 
netic field varies: 0° (top left), 30° (top right), 45° (bottom left), 
60° (bottom right). As the inclination increases, the movement 
of the transition region becomes more regular and also gains in 
amplitude. The horizontal axis is distance along the field rather 
than height, in order to show the true length of the fibrils. The 
left boundary is at a different coordinate in each panel, so that the 
unperturbed transition region appears in the same position, but all 
panels are at the same scale. 

even stronger piston driving we have produced decelera- 
tions that are slightly greater than gravity. Such decel- 
erations have also re cently been reported in observation s 
of quiet sun mottles (|Rouppe van der Voort et aT1 l2007). 

One difference betw een our results and those of 
iDe Pontieu et all (|2QQ7l ) is that they find many more fast 
fibrils than we do. Significant numbers have maximum 
velocities up to 25-30 km s _1 , whereas we only find such 
values for short-period waves propagating vertically. In 
fact, there are few fibrils with greater than 20 km s _1 
maximum velocity in ou r simulations. It is also worth 
noting that the data of IDe Pontieu et all are not cor- 
rected for projection effects. We will touch on a possible 
explanation for this discrepancy later in this section. 

We also find fibrils with lower maximum velocity than 
the lowest found in the observations. This is probably be- 
cause the fibrils are quite easy to find and isolate in our 
simulation data, while superposition effects and back- 
ground noise make small disturbances hard to detect in 
observations. These fibrils show a clear correspondence 
with the wave movement in the simulations, so they are 
not an effect of background noise there. 

We now turn our attention to the question of whether 
the inclination of the field helps lo ng-period wave s prop - 
agat e and reach c orona l heights. iMichalitsanosI (|1973l ) 
and iBel fc Lerovl (|l977l ) made early investigations into 
this phenomenon. The idea is that if the plasma is con- 
fined to moving along the magnetic field, it will also be 
subjected to a reduced effective gravity g cos 0, where 6 is 
the inclination of the magnetic field. With the acoustic 
cutoff period given as 
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where 7 is the ratio of specific heats and c s the 
sound speed, a reduced effective gravity also leads to a 
higher cutoff period, thus allowin g longer-period waves 
to propagate. More r ecentl y, ISuematsTu ([T990) and 
iDe Pontieu et all (|2QQ4l . l2QQ5h have proposed that this 
mechanism may let solar p-modes, with a period of 
around 5 minutes, propagate through the chromosphere 
and provide the driving force for spicules and fibrils. 
These waves are usually evanescent, having periods 
above the acoustic cutoff. In our model atmosphere, the 
cutoff period at the lower boundary is 213 s at 0° in- 
clination, 245 s at 30°, 301 s at 45°, and 425 s at 60°. 
The real solar chromosphere should have values close to 
these. 

Note that the reduced effective gravity hypothesis only 
works if the plasma is in fact confined to moving along 
the magnetic field, that is, in \ow-f3 plasma. In our model 
this applies everywhere, but on the actual sun, where the 
main source of wave energy is the convection region and 
photosphere, we would only expect long-period waves to 
be able to propagate upwards from the photosphere with 
significant energy in areas where the magnetic field is 
strong enough to dominate over thermal pressure forces 
even at photospheric depths. This applies mainly in the 
vicinity of sunspots, active regions and network flux con- 
centrations. 

In Fig. [TT] we have plotted the movement of the transi- 
tion region with time for four different inclinations, with 
the other parameters — piston period and amplitude — 
constant. The period is 300 s, and waves with this period 
ought to be evanescent in a vertical field, being above the 
cutoff period of 213 s. 

With a vertical field (top left), there is some tunneling 
of wave energy, but the movement is slight and some- 
what irregular. As the inclination of the field increases 
to 30° (top right), 45° (bottom left) and 60° (bottom 
right), there is a clear progression towards more regular 
parabolic movement, and the amplitude also increases 
m arkedly. T hese results clearly support the hyp othesis 
of lSuematsul ([l990h and lDe Pontieu et all (I20M l2005h . 
They also support the claim by IDe Pontieu et al.l (|2007l ) 
that the regional differences in behaviour they observe 
are at least partially caused by the inclination of the 
magnetic field, allowing longer period waves (with cor- 
responding lower fibril decelerations) to propagate in re- 
gions of inclined field. 

In Table [TJ we have listed the averaged maximum ve- 
locities of all fibrils at given magnetic field inclinations, 
driver amplitudes and periods. Unsurprisingly, at 0° in- 
clination, only 180 s waves produce proper shocks, al- 
though it is possible to produce regular movement of the 
transition region with 240 s waves if the driving piston 
is strong enough. There is of course a dependence of the 
maximum velocity on the driver amplitude, but it is not 
linear. The maximum velocity at 180 s period increases 
by only about 10% when the driver amplitude increases 
from 0.8 to 1.1 km s _1 , a pattern that is repeated for 
several other combinations of period and inclination. It 
also increases by only a factor of 2.4 when the amplitude 
increases by a factor of 5.5, from 0.2 to 1.1 km s _1 . This 
indicates that the waves are reaching a plateau where 
larger initial amplitudes just lead to increased dissipa- 
tion. 

The maximum velocity of 180 s waves decreases when 



TABLE 1 

Maximum velocities (km s _1 ) 



0° inclination 






Driver amplitude (km s 1 ) 




Period (s) 


0.2 


0.5 


0.8 


1.1 


180 


10.6 


18.7 


22.2 


24.0 


240 






6.2 


8.4 


300 










360 










30° inclination 


180 


9.9 


17.1 


18.7 


19.2 


240 




6.1 


9.7 


12.9 


300 










360 










45° inclination 


180 


5.0 


11.1 


13.5 


14.8 


240 


6.6 


15.3 


18.7 


20.0 


300 






7.1 


9.9 


360 










60° inclination 


180 










240 




8.0 


11.2 


13.2 


300 




10.2 


14.0 


15.5 


360 




8.8 


14.1 


17.1 



Note. — The values are averaged maximum velocities for all 
fibrils with that driver period and amplitude, corrected for projec- 
tion effects. At the missing points, the movement of the transition 
region is irregular or has very low amplitude (j 5 km s — 1 ). 



the inclination increases. Since these waves are well be- 
low the acoustic cutoff, this can not be because they are 
becoming evanescent; instead, it is likely because of the 
increased effective travel distance when the field is in- 
clined. At 60° inclination, the effective distance to the 
transition region is twice as long. This leads to increased 
dissipative losses as the wave travels through the chromo- 
sphere. The effect is more pronounced at 180 s because 
short-period waves steepen into shocks sooner and there- 
fore are more vulnerable to dissipation. 

At the other periods, the maximum velocities increase 
as the acoustic cutoff period increases, becoming notice- 
able when the wave period is similar to the cutoff and 
significant when it is well below. At 60° inclination, the 
240 s waves drop in amplitude again, while those with 
longer periods are at their strongest. 

It should be noted that field inclination is not the only 
possible way to increase the cutoff period. We see in 
equation [2] that we can also modify the sound speed, 
which is proportional to the square root of the tempera- 
ture. Hence, in a locally hotter medium, we can get easier 
propagation of long-period waves. This mechanism may 
be unlikely to be enough on its own, as you need a tem- 
perature increase by a factor 2 to get the same cutoff 
increase as a field inclination of 45°, and a factor 4 to 
match 60°; however, the combination of the two mech- 
anisms could be very effective. For example, in the real 
solar atmosphere, the field tends to be more vertical at 
lower heights but becomes more inclined higher up due 
to the natural expansion of flux tubes. It is then possible 
that a local temperature enhancement could increase the 
cutoff period in the lower layers while the field inclination 
takes over as the waves propagate upwards. Preliminary 
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experiments using a simple model with vertical field up 
to about 1000 km and inclined fields above that indicate 
that such a configuration can also let long-period waves 
propagate. 

It could also be a possible reason why, in our simula- 
tions, we see fewer fibrils with very high maximum ve- 
locities than lDe Pontieu et al.l (|2QQ7h . In our model, the 
field inclination is constant, leading to a large increase in 
effective travel distance at higher inclinations. In a more 
realistic model with the field inclination increasing grad- 
ually, it is possible that the waves travel more vertically 
through lower parts of the atmosphere before the field 
inclination increases, leading to shorter travel distances 
and less dissipation. 

4. SUMMARY 

Our simulations have shown that, even with a simple 
ID model, we are able to reproduce the main observed 
properties of dynamic fibrils. We get parabolic shapes 
and reproduce the range of deceleratio ns and roughly the 
range of maximum velocities found by iDe Pontieu et all 
(2007). We find that the slope of the correlation between 
maximum velocity and deceleration varies with the driver 



period, with lower decelerations at longer periods. The 
distribution of decelerations is incompatible with a bal- 
listic model but fits very well with a shock wave decel- 
eration model. Furthermore, we have shown that long- 
period waves can propagate and reach the transition re- 
gio n if they travel in a st rong i nclined field, as sug gested 
by iMichalitsanosI (| 1973( 1 and iBel fc Lerovl tlflZa), and 
that they can drive fibrils there as su ggeste d bv lSuematsul 
(j!990h and lDe Pontieu et all (|2QQ4J . l2QQ5h . The different 
slopes of the correlation and the leakage of long-period 
waves can account for the regional differences in the ob- 
servations. Our results give strong support to the theory 
that jets such as dynamic fibrils and quiet sun mottles 
are driven by slow mode magnetoacoustic shocks in the 
chromosphere. 
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